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ABSTRACT 

A set 5 of n points in general position in R'' defines the unique Voronoi 
diagram of S. Its dual tessellation is the Delaunay triangulation (DT) of S. In 
this paper we consider the parabolic functional on the set of triangulations of S 
and prove that it attains its minimum at DT in all dimensions. The Delaunay 
triangulation of S is corresponding to a vertex of the secondary polytope of 
S. We proposed an algorithm for DT's construction, where the parabolic func- 
tional and the secondary polytope are used. Finally, we considered a discrete 
analog of the Dirichlct functional. DT is optimal for this functional only in two 
dimensions. 

Keywords. Voronoi diagram, Delaunay triangulation, regular triangula- 
tion, secondary polytope, flip, incremental algorithm. 

1 Introduction 

Some of the most well-known names in Computational Geometry are those of 
two prominent Russian mathematicians: Georgy F. Voronoi and Boris N. De- 
launay. Their considerable contribution to the Number Theory and Geometry 
is well known to the specialists in these fields. Surprisingly, their names (their 
works remained unreaded and later re-discovered) became the most popular not 
among "pure" mathematicians, but among the researchers who used geometric 
applications. Such terms as "Voronoi diagram" and "Delaunay triangulation" 
are very important not only for Computational Geometry, but also for Geomet- 
ric Modeling, Image Processing, CAD, Geographic Information System etc. 

The Voronoi diagram is generated by a set of n points S = {xi, ...,a;„} in 
R"*. The Voronoi diagram is the partition of the R"^ into n convex cells, the 
Voronoi cells Vi, where each Vi contains all points of the R'' closer to Xi than 

*The research was partially supported by the RFFI grant 97-01-00174 and by the CTIAC 
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to any other point: 

Vi = {a;|Vj 7^ i,d{x,Xi) < d{x,Xj)}, 

where d{x, y) is the Eudidean distance between x and y. 

This concept has independently appeared in various fields of science. The 
earliest significant use of Voronoi diagrams seems to have occurred in the work 
of Gauss, Dirichlet and Voronoi in their investigations on the reducibility of 
positive definite quadratic forms. Different names particular to the respective 
field have been used, such as medial axis transform in biology and physiology, 
Wigner-Seitz zones in chemistry and physics, domain of action in crystallogra- 
phy, and Thiessen polygons in meteorology and geography [1]. 

For generic set of n points S in R'' the straight-line dual of the Voronoi 
diagram is triangulation of S, called the Delaunay triangulation and denoted 
by DT{S). The DT{S) is triangulation of the convex hull of S in R"^ and set of 
vertices of DT{S) is S. 

Voronoi [21] was the first to consider the dual structure of the Voronoi dia- 
gram, where any two points of S arc connected whose regions have a boundary 
in common. Later Delaunay [6] obtained the same by defining that two points 
of S are connected if and only if they lie on a sphere whose interior contains 
no points of S. After him, the dual of the Voronoi diagram has been denoted 
Delaunay tessellation or Delaunay triangulation. 

Voronoi diagrams and Delaunay triangulations are used in numerous appli- 
cations. It is widely used in plane and 3D case. A natural question may arise: 
why those structures are better than the others. Usually the advantages of pla- 
nar Delaunay triangulation are rationalized by the max-min angle criterion and 
other properties [1,2,3,6,16,17,19]. 

The max-min angle criterion requires that the diagonal of every convex 
quadrilateral occurring in the triangulation "should be well chosen" [19], in 
the sense that replacement of the chosen diagonal by the alternative one must 
not increase the minimum of the six angles in the two triangles making up the 
quadrilateral. Thus the Delaunay triangulation of a planar point set maximizes 
the minimum angle in any triangle. More specifically, the sequence of triangle 
angles, sorted from sharpest to least sharp, is lexicographically maximized over 
all such sequences constructed from triangulation of S. 

In the papers [13 15] we defined several functionals on the set of all trian- 
gulations of S in R-^ attaining global minimum on the Delaunay triangulation. 

The "mean radius" functional is the mean of circumradii of triangles for pla- 
nar triangulations. Let t be a triangulation of S in the plane. Assume that each 
triangle Aj of this triangulation is related to the radius Ri of its circumcircle. 
Thus every triangulation t is related to the set {i?Ai, -^Afc} of circumradii 
of triangles G t. The numbers of triangles for any two triangulations of S 
are equal, so it is possible to compare sets of radii for different triangulations. 
In particular, it is possible to compare sums of radii: J2^^i power sums: 
^ i?^. , a > 0. It seems that triangulation having minimal sum of radii is "bet- 
ter", because all its triangles in "average" are nearer to the regular triangles. 
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The functional R(t,a) = ^i?^,,a > attains its minimum iff t is Delaunay 
triangulation [15]. 

The harmonic index of triangulation has its origin in the theory of the so 
called "harmonic maps". For polygon P its harmonic index 



where oi, . . . , a„ are the lengths of sides of P and S{P) is its area. This index 
is the same for similar polygons. It is easy to prove that harmonic index for 
triangles achieves its minimum iff triangle is equiangular. For planar triangula- 
tion t of a set 5 let denote by hrm{t) (harmonic index of triangulation t) the 
sum of hrm of its triangles: 



. Harmonic index hrm{i) of triangulation t oi S achieves its minimum iff t is 
the Delaunay triangulation of S. The harmonic functional for triangle attains 
its minimum if this triangle is equiangular. Usually, a triangulation is regarded 
as "good" for different purposes if its triangles are nearly equiangular. The 
harmonic index of triangulation t achieves its global minimum if i is a part of 
a regular triangular lattice. In some sense, this result shows that the Delaunay 
triangulation is as close as possible to equiangular triangulation. 

One of the most popular algorithm for constructing planar Voronoi diagram 
- Delaunay triangulation is so called flipping algorithm [1,8,10]. (If t is planar 
triangulation of 5* and AC is an internal edge of t, the two triangles ACB 
and ACD of t incident with AC in triangulation. If quadrilateral ABCD is 
convex a new triangulation of S may be obtain by removing the edge AC and 
inserting the edge BD. This operation is called the flip or edge-flip). The 
proof that sequence of flips does not cycle in the algorithm easily follows from 
consideration of the mean-radius or harmonic functional. Indeed, using flipping 
algorithm after each flip the functional decreases until Delaunay triangulation 
is reached. 

Among the various proposed methods for constructing Voronoi diagrams in 

R'', incremental insertion of points is most intuitive and easy implement. Joe 
[9], Rajan [16], and Edelsbrunner and Shah [8] generalized of planar incremental 
and flipping algorithms for higher dimensions. 

This paper includes algorithm for construction of the Voronoi diagram. First, 
in section 1 we consider a so called parabolic functional on a set of all triangula- 
tions of S and prove that this functional attains its minimum on the Delaunay 
triangulation of S in all dimensions. The secondary polytope is the original one 
due Gel'fand, Kapranov, and Zelevinsky [4,5]. They introduced the secondary 
polj^ope Q oi S d R'' that is a convex polytope in R",n = ji?!, and where the 
vertices of Q are in one-to-one correspondence with the regular triangulations of 
S. In section 2 we show that parabolic functional is a linear function on R", and 
Delaunay triangulation of 5 is a vertex of Q that gives minimum for parabolic 
function. The incremental algorithm for construction of the Voronoi diagram i.e. 



hrm{P) = Y,aVS{P) 
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Delaunay triangulation is given in section 3. The main idea of this algorithm is 
to construct a sequence of a regular triangulations (vertices of Q) that decrease 
parabolic functional, and the last triangulation in this sequence is DT. Finally, 
in section 4 wc consider a discrete analoque of the Dirichlct functional on a set 
of all triangulations of S. For d = 2 this functional achieves its minimum for DT. 
We consider case n = d + 2 and prove that optimal triangulation depends on 
points S configuration only. If d > 2 then Delaunay triangulation is not optimal 
for Dirichlet functional. Thus the problem to finding "good" triangulations for 
this functional in higher dimensions is opened and more detailed consideration 
is necessary. 

2 Optimality of Delaunay triangulations for the 
parabolic functional 

Throughout this paper S = {xi,...,a;„} denotes a set of n points in general 
position in R''. A triangulation of is a triangulation of the poly tope P = 
CH{S) (convex hull of S) with vertices in S. 

Let f be a triangulation of the set S in R'', denotes the i-th d-simplex of 
t and Xij € R'*, j = 0, 1, . . . , d are its vertices . Let 

Vr{t)=Y,{xlo + -+xl>ol{Ai), 

i 

where vol{Ai) is volume (area for d = 2) of the simplex A^. We call functional 
Vr parabolic (or Voronoi). 

The parabolic functional induces an order on triangulations of the set S by 
the rule: ti > t^ iff Vr{ti) > Vrit^)- The value of Vr depends on the choice of 
the origin. If we move the origin to xq then this order does not change i.e. 

Vr{t{) > Vr{t2) iff Vr{ti,xo) > Vr{t2,xo). 

The main result for Voronoi functional Vr is the following: 

Theorem 1. The parabolic (Voronoi) functional Vr{t) achieves its minimum 
if and only if t is the Delaunay triangulation 

A simple proof of these theorems follows from paroboloid construction of the 
DT found by Voronoi [21] and rediscovered only in 1979 for a sphere (K.Brown), 
and later also for a paraboloid. 

Let us consider an arbitrary triangulation t of the set S and "lift" it onto the 
paraboloid in R''+^, i.e. let us build a polyhedral surface in R''^^ connecting 
corresponding vertices on the paraboloid. Note that the functional Vr up to a 
constant equals to the volume of the solid body below this surface. Thus, the 
minimum of Vr is attained on the Delaunay triangulation. 

Let us consider another functional for triangulations: 

C2(0 = ^||c(Ai)|p^;o;(A,), 
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where c(Ai) is the center (barycenter) of the Aj, c(Aj) = Xij/{d + 1) and 
1 1 • 1 1 is Euchdean norm. 

By direct calculation (it is sufficient to check the formula on a simplex) it is 
possible to prove that 

{d+lfC2{t) + Vr{t) = id+l)id + 2) J \\x\\^dx, 

CH{S) 

where CH{S) is convex hull of set S in R''. 

From Theorem 1 and this formula directly follows that: 

Theorem 2. The functional C2(t) on triangulations of the set S achieves its 
maximum if and only if t is the Delaunay triangulation. 



3 The secondary polytope, regular triangulations, 
and flips 

The following analytic description of the secondary polytope is the original one 
due Gel'fand, Kapranov, and Zelevinsky [4,5]. We include it here for complete- 
ness. They introduced the secondary polytope Q = X^(<S') of an afSne point 
configuration S, where the vertices of ^-rc in one-to-one correspondence 

with the regular triangulations of the "primary polytope" P = CH{S) - convex 
hull of S. 

Let S = {xi} be a set of n points in R'' and t is a triangulation of S. We 
are correspond to triangulation t vector p{t) = {pi,p2, ■ ■ ■ ,Pn) in R-", where 
Pi = Y1 vol{Aij, Aij denotes the j-th d-simplex of t that incident to vertex Xi, 
and vol{A) denotes volume of the d-simplex i.e. Pi is volume of star of i-th 
vertex of triangulation t. 

We have d+1 equations: 

= {d+l)vol{CH{S)), ^XiPi = {d+l)xcVol{CH{S)), 

where Xc is center (center of mass) of CH{S). Since the right hand sides of these 
equations not depend of triangulation that in fact p{t) G R"~'^~^. Let 

Q = convex hull{p(f) :t is a triangulation of S}, 

i.e. Q is a convex hull of the set of images in R" all triangulations of S. 
The convex polytope Q called secondary polytope and denoted by ^{S). The 
dimension of this polytope is n — d — 1. 

Even for a simple configuration of S the secondary polytope is not simple. 
When S in the plane consist of n vertices of a convex n-gon, then secondary 
polytope called associahedron [12] . Associahedron is a simple (n-3)-dimensional 
polytope with n{n — 3)/2 facets. For example, associahedron of a pentagon 
is a pentagon, and associahedron of a hexagon is a simple 3-polytope with 14 
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vertices, 21 edges, and 9 facets. An interesting application of the associahe- 
dron to the theoretical Computer Science has been given by Sleator, Tarjan, 
and Thurston [20]. These authors derive a tight upper bound for the rotation 
distance between binary trees with n nodes by proving that the diameter of the 
associahedron equals 2n — 10, for large n. 

A triangulation f of 5 is said to be regular if there exists a function on P that 
is piecewisc linear and strictly convex with respect to t. (A convex piecewise 
linear function over a triangulation t is said to be strictly convex if it is given 
by a different linear function on each maximal cell of t). 

Regular triangulations are really just the duals of power diagrams, and Edels- 
brunner call them weighted Delaunay triangulations [8] . 

There are several equivalent ways to define the notion of a regular triangu- 
lation of S. For example, Gel'fand, Kapranov, and Zelevinski [5] call regular a 
following triangulation t: 

Choose numbers yi,. . . ,yn and let W = {(a;i, yi), . . . , {xn, Vn)} C R'^+^. If 
A = {{xi^ , j/i J, . . . , (.Ti^^i , yid+i)} is a facet of W in the lower hulloi W (i.e., the 
last component of the outward normal of its supporting hyperplane is negative) 
then A = {a;^^ , . . . , Xi^^^ } is a d-face of the triangulation t. 

Let yi = Then we get the Delaunay triangulation of S. Therefore, 

from this definition the Delaunay triangulation is regular. 

Consider a set 5 of rf + 2 points in R"^. From Gale diagram [7] follows (see 
also Schlcgcl [18], Lawson [11]) that there are exactly two ways to triangulate 
S. Indeed, the two ways correspond to the two sides (lower and upper) of the 
{d + l)-simplex that is convex hull of corresponding lifted points in R'*+^. A 
flip is the operation that substitutes one triangulation of S for the other [8] . 

There are three types of flips in two dimensions, and we denote a flip by the 
number of triangles before and after the flip. So the flips in two dimension are 
of type '1 to 3', '2 to 2', and '3 to 1'. The first type introduces a new point, and 
the last type removes a point. The flips in three dimensions can be classified as 
'4 to 1', '3 to 2', '2 to 3', and '1 to 4'. 

Now we can give the main results of this section. First theorem belongs to 
Gel'fand, Kapranov, and Zelevinsky. They used another terminology. 

Theorem 3. The vertices of the secondary polytope Q = ^{S) are in one-to- 
one correspondence with the regular triangulations of S, and the edges of Q are 
corresponding to flips. 

From this theorem and Section 1 follow: 

Theorem 4. There is a sequence of flips that connect any regular triangulation 
of S with the Delaunay triangulation of S, and after each flip decreases parabolic 

functional. 

Proof. It is easy to see that for any triangulation t a parabolic functional Vr{t) = 
J2 IkilPPi) where pi = J2'^ol{Aij) as above, pi is i-th coordinate of R" and 
therefore Vi — functional is a linear function on R". Secondary polytope Q = 
^(5) is a convex polytope in R", and Delaunay triangulation of S {DT{S)) is a 
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vertex of Q that gives minimum for function Vr. It is clear how to find sequence 
of neighboring vertices of Q decreasing parabohc functional and connected any 
regular triangulation (vertex of Q) with vertex corresponding to the Delaunay 
triangulation. □ 



4 Incremental construction of the Voronoi dia- 
gram. 

A natural idea is to construct the Voronoi diagram by incremental insertion, i.e. 
to obtain V{S) from V{S)\{x} by inserting the point x. The insertion process 
is, maybe, better described, and implemented in the dual environment, for the 
Delaunay triangulation: construct DTi = DT{{xi, X2, Xi^i, Xi}) by inserting 
the point Xi into £)Tj_i. The advantage over a direct construction of V{S) is 
that Voronoi vertices that appear in intermediate diagrams but not in the final 
one need not be constructed and stored. 

Several algorithms proposed for Delaunay triangulation are based on the 
notion of a local transformation henceforth referred to as a flip. Historically 
the first such algorithm is due to Lawson [10]. Given a finite point set in the 
plane, the algorithm first construct an arbitrary triangulation of the set. This 
triangulation is then gradually altered through a sequence of edge-flips until the 
Delaunay triangulation is obtained. The generalization of this method to R'', 
d > 2 has difficulties, and it is incorrect if the fiips are applied to an arbitrary 
initial triangulation [8]. Joe [9] shows that if a single point, Xi, is added to the 
Delaunay triangulation DTi^i in R"^ then many different sequences of flips will 
succeed in constructing the Delaunay triangulation DTi. This can be used as the 
basis of an incremental algorithm. Rajan [16] considers Delaunay triangulation 
in arbitrary dimensions, R**, and argues that a single point can always be added 
by a sequence of flips. However, he needs a priority queue to find the appropriate 
sequence, which takes logarithmic time per flip. Edelsbrunner and Shah [8] 
using "weighted points" method unifies and extends the algorithmic rcsiilts of 
Joe [9] and Rajan [16]. In particular, they show that many different sequences 
of flips can be used to add a single point to a regular triangulation in R"^. This 
eliminates the need for a priority queue that sorts the flips. This section wc show 
how to applied parabolic functional and secondary polytope for construction of 
the Voronoi diagram incrementally. 

The algorithm constructs the Delaunay triangulation of a give set S = 
{xi,X2, ■■■,Xn} in R'' incrementally. It is convenient to first construct an ar- 
tificial rf-simplex 5*0 = {xd, a;_rf+i, xq}, so that S is contained in it. The 
d + I artificial points can be conveniently chosen at infinity, so that choice 
of points guarantees that DT{S) is a subcomplex of DT{S[j So) [8]. In fact, 
DT{S) consist of all simplices of DT{S{J So) that are not incident to any point 
of S-Q. 

Let i be a triangulation of S. Call t locally non-optimal triangulation if 
there is a convex subcomplex a & t that after the flip a a new triangulation is 
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regular and decreases parabolic functional (sec section 1.) i.e. Vr{t) > Vr{t'), 
where t' is the new triangulation. In accordance with this definition we call 
t L.O.T. (locally optimal triangulation) if t is not locally non-optimal. We 
denote a locally optimal triangulation of S as LOT{S). We will show later that 

LOT (Si) = DTi for alH = 0, 1, n. 

THE INCREMENTAL ALGORITHM FOR CONSTRUCTION OF LOT{S) 

1 Construct LOT{So); 

2 for i := 1 to n do 

3 locate the d-simplex s in LOT{Si-i) that contains Xi 

4 if iOT(5i_i \J{xi} is not L.O.T. then 

5 flip LOT{S^-i[}{xi}- 

6 while there exist non-optimal d + 2-subcomlex do 

7 find a non-optimal d + 2-subcomplex cr; 

8 flip a 

9 endwhile 

10 endif 

11 endfor 

This algorithm could fail for two reason. First, it could be that the while 
loop does not terminate, because it cycles in an infinite loop of flips. Second, 
if the algorithm would stop before reaching the DT{S). We show this cannot 
happen. 

The proof that sequence of fiips in the algorithm does not cycle easily follows 
from the consideration of the parabolic functional Vr{t). Indeed, this functional 
decreases after each flip in the algorithm. 

For each i in algorithm LOT{Si) is locally optimal triangulation of Si. It 
means, that LOT {Si) is regular triangulation, and there is no a locally non- 
optimal subcomplex a. Then LOT(Si) is a vertex of a secondary polytope 
Qi = J2{Si)- DTi is vertex of Qi also, and therefore from Theorem 4 (see sec- 
tion 2.) follows that there is sequence of flips that connect LOT{Si) and DTi. 
From theorem 1 follows that DTi gives minimum for parabolic functional. From 
other side parabolic functional cannot be decrease for LOT{Si). Consequently 
LOT{Si) is DTi, and therefore 

Theorem 5. Constructed in the algorithm a locally optimal triangulation of S 
(LOT{S) ) is the Delaunay triangulation of S. 

5 The Dirichlet functional on triangulations 

Let S = {xi} be a set of n points in R'^', each associated with a real number yi. 
Denote by Y the set of these numbers, i.e. Y = {yi, ...,yn). There are a lot of 
different problems in Geography, Geology, Topography, CAD/CAM etc., where 
we need to construct a surface in R''"^^ corresponding to this dataset. The main 
problem is the following: to find a function y = f{x), such that f{xi) = yi. One 
of the oldest and the most famous methods is modeling by triangulation. If we 
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have some triangulation of S then for a set of data Y there is only one method to 
construct a piecewise hnear function (polyhedral surface) on this triangulation. 
Usually Delaunay triangulation is used for this purpose. 

One of the minimum criterion is a discrete analogue of the Dirichlet func- 
tional: / \\grad f{x)\\'^dx. For interval {d = 1), a spline of deg = 2A; — 1 is a 
function y = f{x) such that /(a;,) = y, and 

6 

j [f''\x)fdx = min. 

a 

For k = \^d> 1 and piecewise linear function / we get 

DF{t,Y)= j \\gradf{x)\fdx = ^ ^""'^l^f^^' ' vol{CH{S)). 

CH{S) ' 

The triangulation of S that minimizing the functional DF can be called the em 
discrete spline triangulation (DST). For the plane DST does not depend upon 
Y and it is DT. Rippa [11] for d = 2 proved that DF{t, Y) achieves its minimum 
iff t is DT. 

For d > 2 the Delaunay triangulation of S could be not optimal for this 
functional. Let 

a;o = (0,0, ••• ,0); Xj = ( 0,-- ,0 ,1, 0,-- ,0) ,i= I,--- ,d; Xd+i = {a,- ■ ■ ,a). 

i—l n—i 

There are exactly two distinct triangulations of 5 = {xo, xi, - ■ ■ , Xd+i)- Denote 
by ti a triangulation consist of two simplices: Ai = {xo,xi,--- ,Xd)\^2 = 
{Xci+i,xi, - ■ ■ ,Xd), and by t2 another one. Then t\ is DT iff a > 1. It is 
easy to show by direct calculation that t\ is DST iff o > Therefore, for 

3^ < a < 1 is DST, but is not DT. 

The proof of Rippa's theorem directly follows from the fact that the DF- 
functional for triangulation f of a quadrilateral is minimum if t is DT. The proof 
also follows from some general result that is given below. 

Let S be a set of d + 2 points xi, ...,Xd+2 in R''. Suppose S admits two 
triangulations ti and t2, and Y = (yi, ...,yd+2) is a set numbers corresponding 
to xi,...,Xd+2 as above. Let B(Y,S) = DF{ti,Y) - DF{t2,Y). Note B{Y,S) 
is a quadratic form depending on Y. 

Theorem 6. The optimal (DST) triangulation of S for n = d + 2 does not 
depend on Y. 

Proof. It is easy to see that for arbitrary set of real numbers (ao, ai, • • • , Od) '■ 
B{Y,S) = B{Y,S), where F = (^i, • ■ • ,y„), and 

d 

Vi = yi + ao + ^ajXij; Xi = {xa,--- ,Xid) G S. 
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Then for d + l-dimensional subspace K C R" that is the linear hull of the 
set of vectors: 

(I,-- - ,1), {Xu,--- ,Xnl),--- ,{xid,--- ,Xnd) 

the quadratic form B{Y, S) is vanished, and B{Y + Z,S) = B{Y, S) if Z e 3?. 
Note that n ~ d + 2 and dimSR = d + 1 therefore B{Y, S) could be not vanish 
only on l-dimensional subspace of R" that is orthogonal to 3?. Thus B{Y, S) = 
const{S)L'^{Y), where L{Y) is some linear form on F £ R". Therefore sign of 
the B{Y,S) = DF{ti,Y)-DF{t2,Y) does not depend on F, and if ii is optimal 
for some Y (i.e. DF{ti,Y) < DF{t2, Y)) that it is optimal for any Y. □ 

We state an open problem concerning DST: 
- Does DST depend on Y, when n> d + 2 ? 

6 Concluding Remarks. 

Voronoi diagram and Delaunay triangulation have a fair number of applications, 
including the generation of grids for point configuration and for surface inter- 
polation. No doubt that in two-dimension these tesselations are the best for 
these purposes. Indeed, the main motivation for studying the problems solved 
in this paper is our intention to implement Voronoi and Dirichlet functionals 
in dimensions beyond R^. The optimal triangulations for these functionals in 
R'', d > 2 could be not the equal. We do not know is optimal triangulation 
for Dirichlet functional is regular? If it is regular then for it construction the 
algorithm in section 3 is suitable. Instead of Fr-functional there have to be 
used DF functional. It would be interesting to study dual tessellation for DST. 
In other word, what is analogue of Voronoi diagram for DST? 
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